Vektorová a rastrová geodata v R, jejich interakce a vykreslování

Autor

Ondřej Ledvinka

Publikováno

9. prosince 2025

Upraveno

15. prosince 2025

Prerekvizity

V následujícím se předpokládá, že již máme založený R projekt ze začátku kurzu. To znamená, že máme možnost využívat relativní cesty k souborům a jiné výhody práce s R projektem. Dále se předpokládá, že máme kvalitní připojení k internetu, jelikož se může stát, že budeme instalovat chybějící balíčky. Internetové připojení je nutností při získávání geodat z internetu, což si také v některých případech předvedeme.

Pokud možno, budeme také pracovat s nejnovějšími verzemi R (v současnosti R version 4.5.2 (2025-10-31 ucrt) – [Not] Part in a Rumble z 2025-10-31) a RStudio (v současnosti RStudio 2025.09.2+418 “Cucumberleaf Sunflower” z 2025-10-20).

Vektorová geodata

Vektorová geodata v R lze ve stručně charakterizovat jako tabulky obohacené o sloupec, ve kterém máme navíc informaci o geometrii, ke které se ostatní atributy (v ostatních sloupcích – polích) váží. Nejvíce budeme pracovat s tzv. simple feature collection (sf), kterou lze ale rozdělit do několika prvků. Samozřejmostí je atributová tabulka nebo jen jedna řádka s atributy (či dokonce jedno pole s jedinou hodnotou), ale zvláštností zde je sloupec s geometrií, tzv. simple feature column (sfc). Tento sloupec se skládá z jedné nebo více simple feature geometry (sfg). Podle těchto tříd pak lze aplikovat funkce nejdůležitějšího R balíčku sf, které tyto třídy akceptují. Kromě toho se můžeme zabývat také vztahy mezi atributy a geometrií (attributes-geometry relationships, agr), ale do toho se pouštět nebudeme1.

Pokud jde o geometrii, rozlišujeme zde tzv. velkou sedmičku – (MULTI)POINTS, (MULTI)LINESTRING, (MULTI)POLYGON a GEOMETRYCOLLECTION. Velice zřídka se setkáme s dalšími, vzácnějšími, typy podporovaných geometrií, přičemž kompletní výčet uvádí kniha Pebesma, Bivand (2023). Tam se také dočteme, že se na typy geometrie lze doptávat2 nebo dokonce můžeme mezi nimi “přepínat”, pokud to logika věci dovoluje (viz funkci st_cast()). S geometrií lze pracovat i tak, že lze přičítat či násobit konstanty (nebo třeba i matice).

Celý tidyverse jsme se vlastně učili proto, že právě jeho přístupy jsou velmi dobře podporovány při práci s třídou sf, a to včetně pipe operátoru, takže nám zde přibývají funkce, které umí pracovat se speciálními typy tabulek s geometrií. Výhodou je, že sf třídu lze funkcí as_tibble() převádět na tabulky a naopak z tabulky můžeme objekt konvertovat na třídu sf pomocí funkcí st_sf() (pokud v podkladové tabulce již existuje sloupec s geometrií, ačkoliv jaksi schovaný) nebo funkcí st_as_sf() (pokud v tabulce neexistuje žádný sloupec s geometrií, ale přesto existuje ve sloupcích informace z níž si geometrii můžeme vytvořit). Funkce st_as_sf() je univerzálnější, podporuje i jiné třídy, než jsou ty tabulkové. Podobně existují funkce vytvářející třídy sfc či sfg.

Ukažme si tedy nějaké příklady ukazující tvorby třídy sfg a tvorbu třídy sf. Dále ukažme, jak lze sf převádět na tabulku a tu opět převádět na sf.

# nejprve musíme načíst potřebné balíčky
# rovnou načteme i balíček pro práci s rastrovými geodaty a další balíček pro případnou kresbu situace
xfun::pkg_attach2("tidyverse",
                  "RCzechia", # důležitý balíček pro vektory sf se načítá automaticky s tímto
                  "terra", # balíček pro rastry
                  "tmap") # pro kreslení map
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.6
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.1     ✔ tibble    3.3.0
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.2.0     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Loading required package: sf

Linking to GEOS 3.13.1, GDAL 3.10.2, PROJ 9.5.1; sf_use_s2() is TRUE

terra 1.8.86


Attaching package: 'terra'


The following object is masked from 'package:tidyr':

    extract
# vytvořme bod z vektoru
# nejpre zadáváme souřadnici související se zeměpisnou délkou
# pak teprve souřadnici související se zeměpisnou šířkou
bod <- st_point(c(15, 50))

# co teď máme za třídu?
bod |> 
  class()
[1] "XY"    "POINT" "sfg"  
# jak objekt vypadá po vytištění do konzole?
bod
POINT (15 50)
# převeďme na třídu sfc, čímž se připravíme tvorbu sf, protože již budeme mít připravený sloupec s geometrií
bod <- bod |> 
  st_sfc(crs = 4326) # zde již můžeme definovat souřadnicový referenční systém (crs)

# číselný kód reprezentuje tzv. EPSG autoritu (zde tedy jde o crs, který známe z GPS, zvaný stručně WGS 84)

# podíváme se, co tu máme teď
bod |> 
  class()
[1] "sfc_POINT" "sfc"      
# tímto se v konzoli objevuje již jakási hlavička
bod
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 15 ymin: 50 xmax: 15 ymax: 50
Geodetic CRS:  WGS 84
POINT (15 50)
# teď ještě převeďme na třídu sf
# nakonec ještě raději přejmenujeme sloupec s geometrií na něco rozumějšího
bod <- bod |> 
  st_sf() |> 
  st_set_geometry("geometry") # tato funkce tady funguje podobně jako set_names(), ale umí toho víc

bod |> 
  class()
[1] "sf"         "data.frame"
bod
Simple feature collection with 1 feature and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 15 ymin: 50 xmax: 15 ymax: 50
Geodetic CRS:  WGS 84
       geometry
1 POINT (15 50)
# protože tohle už je objekt třídy sf, lze přidávat atributy, a to např. funkcí mutate()
bod <- bod |> 
  mutate(nm = "Kouřim")

bod
Simple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 15 ymin: 50 xmax: 15 ymax: 50
Geodetic CRS:  WGS 84
       geometry     nm
1 POINT (15 50) Kouřim
# změňme vykreslovací mód balíčku tmap na interaktivní (tj. na "view")
tmap_mode("view")
ℹ tmap modes "plot" - "view"
ℹ toggle with `tmap::ttm()`
# a nakresleme bod
# funkce balíčku tmap staví na podobné logice jako funkce balíčku ggplot2
# takže pro vrstvení používáme operátor +
tm_shape(bod) + # nejprve definujeme novou vrstvu
  tm_symbols() # pak způsob kreslení
# demonstrujme konverzi třídy sf na třídu tibble
bod <- bod |> 
  as_tibble()

# ztratila se hlavička, ale geometrický sloupec je stále přítomen
bod
# A tibble: 1 × 2
     geometry nm    
  <POINT [°]> <chr> 
1     (15 50) Kouřim
# proto můžeme aplikovat funkci st_sf() při konverzi zpět na třídu sf
bod <- bod |> 
  st_sf()

# teď na tom dokonce budeme ještě lépe, protože tabulka je třídy tibble
bod
Simple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 15 ymin: 50 xmax: 15 ymax: 50
Geodetic CRS:  WGS 84
# A tibble: 1 × 2
     geometry nm    
  <POINT [°]> <chr> 
1     (15 50) Kouřim
# třídy sf se ale můžeme zbavit i jinak, prostřednictvím funkce st_drop_geometry()
# objekt nepřepisujeme, jen demonstrujeme
bod |> 
  st_drop_geometry()
# A tibble: 1 × 1
  nm    
* <chr> 
1 Kouřim

Lze samozřejmě konstruovat složitější geometrické útvary, kde namísto vektorů využíváme matic. Existují i jiné balíčky, usnadňující konstrukci (viz např. sfheaders nebo mapedit, který akceptuje i trvorbu i editaci interaktivním způsobem). Často ale budeme mít k dispozici data se souřadnicemi bodů nebo dokonce již hotová vektorová geodata, která bude postačovat jen načíst. Ukažme si tedy tvrobu bodové vektorové vrstvy ze souřadnic a následně některé možnosti načítání hotových vektorových godat.

# mějme např naši metadatovou tabulku se stanicemi povrchových vod měřícími množství
# tabulku již načíst z JSON souboru umíme
url <- "https://opendata.chmi.cz/hydrology/historical/metadata/meta1.json"

meta <- jsonlite::fromJSON(url)

meta <- meta$data$data$values |> 
  as.data.frame() |> 
  as_tibble() |> 
  set_names(meta$data$data$header |> 
              str_split(",") |> 
              unlist()) |> 
  janitor::clean_names()

# souřadnice jsou teď uvedeny jako textové řetězce
# budeme potřebovat aplikovat funkci st_as_sf(), která ale chce souřadnice jako čísla
meta <- meta |> 
  mutate(across(starts_with("geogr"),
                as.numeric)) |> 
  st_as_sf(coords = c("geogr2", "geogr1"), # tady pozor na záměnu pořadí souřadnic
           crs = 4326)

# protože se chceme zaměřit na geometrický sloupec, raději vybereme jen některé, abychom geometrii viděli
# geometrický sloupec zůstává, i když jsme ho funkcí select() nevybrali
# jedná se o tzv. 'sticky geometry', které se zbavujeme jinak (a my už víme, co pomůže)
meta |> 
  select(dbc:stream_name)
Simple feature collection with 831 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 12.37578 ymin: 48.6255 xmax: 18.8926 ymax: 50.97818
Geodetic CRS:  WGS 84
# A tibble: 831 × 4
   dbc    station_name   stream_name                geometry
   <chr>  <chr>          <chr>                   <POINT [°]>
 1 206200 Slapany        Odrava          (12.37578 50.02859)
 2 293500 Rychvald       Stružka         (18.35676 49.86411)
 3 176100 VD Hracholusky Mže               (13.17649 49.789)
 4 162200 Červená Řečice Trnava          (15.17368 49.52596)
 5 383000 Horní Bečva    Rožnovská Bečva (18.31154 49.42246)
 6 206500 VD Jesenice    Odrava          (12.47487 50.09001)
 7 023500 Orlické Záhoří Divoká Orlice   (16.48049 50.27386)
 8 059000 Nemošice       Chrudimka       (15.78995 50.01647)
 9 070000 Nový Bydžov    Cidlina          (15.4998 50.23809)
10 082000 Plaňany        Výrovka           (15.01614 50.049)
# ℹ 821 more rows
# klidně i vykreslíme
tm_shape(meta) + 
  tm_symbols()
# stáhněme si k tomu např. polygony krajů Česka
kraje <- kraje() # funkce pochází z balíčku RCzechia
RCzechia: downloading remote dataset.
kraje
Simple feature collection with 14 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 12.09057 ymin: 48.55181 xmax: 18.85926 ymax: 51.0557
Geodetic CRS:  WGS 84
# A tibble: 14 × 4
# Groups:   KOD_KRAJ, KOD_CZNUTS3 [14]
   KOD_KRAJ KOD_CZNUTS3 NAZ_CZNUTS3                                     geometry
   <chr>    <chr>       <chr>                                      <POLYGON [°]>
 1 3018     CZ010       Hlavní město Praha   ((14.49219 50.17235, 14.47952 50.1…
 2 3026     CZ020       Středočeský kraj     ((15.14068 50.52093, 15.14296 50.5…
 3 3034     CZ031       Jihočeský kraj       ((14.84258 49.58029, 14.84104 49.5…
 4 3042     CZ032       Plzeňský kraj        ((12.56222 49.6091, 12.56244 49.60…
 5 3051     CZ041       Karlovarský kraj     ((12.95711 49.97298, 12.95698 49.9…
 6 3069     CZ042       Ústecký kraj         ((13.71237 50.7199, 13.71227 50.71…
 7 3077     CZ051       Liberecký kraj       ((15.0068 50.90721, 15.00565 50.90…
 8 3085     CZ052       Královéhradecký kraj ((16.00303 50.65935, 16.00195 50.6…
 9 3093     CZ053       Pardubický kraj      ((16.54661 50.14725, 16.5463 50.14…
10 3107     CZ063       Kraj Vysočina        ((16.21716 49.35523, 16.21751 49.3…
11 3115     CZ064       Jihomoravský kraj    ((17.1357 49.21512, 17.13819 49.21…
12 3123     CZ071       Olomoucký kraj       ((17.04949 49.36472, 17.04903 49.3…
13 3131     CZ072       Zlínský kraj         ((17.77804 49.4432, 17.7785 49.442…
14 3140     CZ080       Moravskoslezský kraj ((17.56338 50.27035, 17.56238 50.2…
# takto můžeme kreslit dvě vrstvy
tm_shape(kraje) + 
  tm_borders() + 
  tm_shape(meta) + # další vrstvu vždy definujeme novým tm_shape()
  tm_symbols()

V balíčku sf také existuje funkce, která vektorová geodata načítá ze souboru. Jmenuje se read_sf() a dokáže načítat z formátů takových, které umožnuje knihovna GDAL. Jelikož tato knihovna umí načítat soubory i bez rozbalení souboru z archivu, nebo rovnou z internetového odkazu, může se povést, že soubory ani nepotřebujeme stahovat, ani robalovat před načtením. Velmi často jsou vektorová geodata také publikována prostřednictvím webových služeb, zde při načítání pomáhají balíčky arcgislayers nebo ows4R. Bohužel není čas toto studovat tak detailně.

Mezi vektorovými vrstvami lze studovat vztahy a na základě těchto vztahů provádět výběry. Nejdůležitější funkcí zde je asi st_intersects(), která je standardně nastavena v dalšíh funkcích, které staví na výběrech, jako je funkce st_join() určená k propojování atributů na bázi prostorových vztahů.

# ukažme si, jak funkci st_intersects() ani nemusíme znát a i tak pomocí ní vybrat body uvnitř nějakého polygonu
# řekněme, že budeme chtít vybrat body jen na území Prahy
# nejprve se na prahu ale budeme muset zaměřit
praha <- kraje |> 
  filter(NAZ_CZNUTS3 == "Hlavní město Praha")

# musíme zajistit stejné crs, můžeme transformovat a dědit z jiné vrstvy
stanice <- meta |> 
  st_transform(st_crs(praha)) # toto je zde ukázáno jen přo případy, kdybychom to potřebovali; zde máme v původních datech již podmínku stejných crs splněnu

stanice_praha <- stanice[praha,]

tm_shape(stanice_praha) + 
  tm_symbols()
# když nastavíme namísto standardně vybrané funkce st_intersects() např. funkci st_disjoint(), dosáhneme pravého opaku při výběru
stanice_nepraha <- stanice[praha, op = st_disjoint]

tm_shape(stanice_nepraha) + 
  tm_symbols()
# podobné výsledky můžeme získat i pomocí funkcí str_crop() a st_intersection()
stanice_praha2 <- stanice |> 
  st_crop(praha) |> # touto funkcí se dostáváme na box daný souřadnicemi v hlavičce ohraničujícího polygonu
  st_intersection(praha) # tato funkce ořezává pak přesněji podle tvaru polygonu
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
tm_shape(stanice_praha2) + 
  tm_symbols()
# výsledek v tomto případě může vypadat stejně, ale zkuste si výsledky porovnat třeba při výběru a ořezávání liniových objektů

# funkce st_join() funguje podobně jako připojování tabulek podle klíčových sloupců, jenom namísto klíčů pracuje na bázi prostorových vztahů
# takto např. můžeme zjistit, kolik vodoměrných stanic máme v jednotlivých krajích
stanice |> 
  st_join(kraje |> 
            select(NAZ_CZNUTS3)) |> 
  st_drop_geometry() |> 
  count(NAZ_CZNUTS3)
# A tibble: 15 × 2
   NAZ_CZNUTS3              n
   <chr>                <int>
 1 Hlavní město Praha      20
 2 Jihomoravský kraj       75
 3 Jihočeský kraj          62
 4 Karlovarský kraj        32
 5 Kraj Vysočina           55
 6 Královéhradecký kraj    61
 7 Liberecký kraj          44
 8 Moravskoslezský kraj   117
 9 Olomoucký kraj          60
10 Pardubický kraj         45
11 Plzeňský kraj           49
12 Středočeský kraj       103
13 Zlínský kraj            57
14 Ústecký kraj            40
15 <NA>                    11
# jak vidíme, existuje i skupina stanic, kterým nebylo možné přiřadit ani jeden kraj

Rastrová geodata

Rastrová geodata jsou většinou moc velká, a proto u nich existují strategie načítání jenom toho, co pro práci v R zrovna potřebujeme. Proto také od začátku kurzu není doporučováno ukládát do a načítat ze souborů typu .RData. Když načteme rastrová geodata, v naprosté většině se Globálním prostředí ukáže, že daný objekt reprezentuje pouze formální třídu. Při tisku do konzole vidíme často jen určitou hlavičku s maximy a minimy hodnot v jednotlivých vrstvách. Více hodnot se z rastrových souborů načítá až tehdy, potřebujeme-li nad nimi provádět výpočty, nebo třeba jen kreslit. Zásadní při zpracování rastrových geodat v R je balíček terra (Hijmans 2025). Pro načítání dat ze souborů používáme funkci rast(). Ukažme si, jak načíst rastrová data ze souborů na internetu, využijme k tomu např. scénářová rastrová geodata vytvořená českými klimatology.

# vybereme pouze teplotu vzduchu a pouze scénář SSP2-4.5
# víme, že na odkaze https://www.perun-klima.cz/scenare/data/SSP245_T_year_asc.zip najdeme ZIP soubor se čtyřmi rastrovými soubory (v ASC formátu)
# názvy souborů můžeme zjistit po stažení a nahlédnutí do ZIP souboru
# využijeme tzv. řetízky a načteme si soubory programaticky
r1 <- rast("/vsizip/vsicurl/https://www.perun-klima.cz/scenare/data/SSP245_T_year_asc.zip/SSP245_T_2021-2040_year.asc")

r2 <- rast("/vsizip/vsicurl/https://www.perun-klima.cz/scenare/data/SSP245_T_year_asc.zip/SSP245_T_2041-2060_year.asc")

r3 <- rast("/vsizip/vsicurl/https://www.perun-klima.cz/scenare/data/SSP245_T_year_asc.zip/SSP245_T_2061-2080_year.asc")

r4 <- rast("/vsizip/vsicurl/https://www.perun-klima.cz/scenare/data/SSP245_T_year_asc.zip/SSP245_T_2081-2100_year.asc")

# funkcí c() můžeme vrstvy spojit, jako bychom tvořili vektor nebo seznam
r <- c(r1, r2, r3, r4)

# podíváme se jak vypadá hlavička objektu při tisku do konzole
r
class       : SpatRaster 
size        : 628, 1068, 4  (nrow, ncol, nlyr)
resolution  : 500, 500  (x, y)
extent      : 263229.1, 797229.1, 5359788, 5673788  (xmin, xmax, ymin, ymax)
coord. ref. :  
sources     : SSP245_T_2021-2040_year.asc  
              SSP245_T_2041-2060_year.asc  
              SSP245_T_2061-2080_year.asc  
              SSP245_T_2081-2100_year.asc  
names       : SSP245_~40_year, SSP245_~60_year, SSP245_~80_year, SSP245_~00_year 
# chybí zde informace o crs, ale my víme, že je to crs s EPSG kódem 32633
# terra potřebuje pro definici crs znát i autoritu
# zde tedy namísto funkce st_crs() existuje funkce crs()
crs(r) <- "EPSG:32633" # akceptovaná jsou i malá písmena

r
class       : SpatRaster 
size        : 628, 1068, 4  (nrow, ncol, nlyr)
resolution  : 500, 500  (x, y)
extent      : 263229.1, 797229.1, 5359788, 5673788  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633) 
sources     : SSP245_T_2021-2040_year.asc  
              SSP245_T_2041-2060_year.asc  
              SSP245_T_2061-2080_year.asc  
              SSP245_T_2081-2100_year.asc  
names       : SSP245_~40_year, SSP245_~60_year, SSP245_~80_year, SSP245_~00_year 
# v balíčku tmap existují i výborné funkce pro kreslení rastrových geodat
# nyní však raději přepneme na kreslení statických map
tmap_mode("plot")
ℹ tmap modes "plot" - "view"
# při výběru vrstvy pro kreslení můžeme použít dvojité hranaté závorky
tm_shape(r[[1]]) + # vybereme průměry za období 2021-2040
  tm_raster(col.scale = tm_scale_continuous(values = "brewer.reds"), # vybíráme lepší barevnou paletu
            col.legend = tm_legend(reverse = T)) # obracíme pořadí barev v legendě
[plot mode] fit legend/component: Some legend items or map compoments do not
fit well, and are therefore rescaled.
ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Funkce rast() toho umí mnohem více a vyplatí se prostudovat její dokumentaci. K ukládání rastrových geodat do souboru naopak slouží funkce writeRaster() nebo writeCDF(), to podle toho, jestli chceme, aby výsledkem byl NetCDF soubor (s příponou .nc; mj. velmi oblíbený v klimatologii)3, či jiný formát, jako je např. GeoTIFF (s příponou .tif). Vypočítejme např. průměrnou teplotu vzduchu pro celé období 2021–2100 a uložme výsledek do GeoTIFF souboru.

# správně bychom měli počítat vážené průměry a vahami by měly počty dnů v jednotlivých obdobích
# váhy již získat umíme, zaměstnejme např. intervaly
vahy <- c(ymd(20210101) %--% ymd(20410101),
          ymd(20410101) %--% ymd(20610101),
          ymd(20610101) %--% ymd(20810101),
          ymd(20810101) %--% ymd(21010101)) / days(1)

vahy
[1] 7305 7305 7305 7304
# nyní můžeme váhy využít pro výpočet váženého průměru rastrových vrstev (počítá se napříč všemi vrstvami, tj. pro každou buňku dostaneme časový průměr)
# k tomuto dobře poslouží funkce terra::weighted.mean()
tavg <- weighted.mean(r, w = vahy) |> 
  round(1) # zaokrouhlujeme, protože nemá smysl uvádět teplotu na více jak jedno desetinné místo

# prohlédněme hlavičku
# je zde skutečně už jenom jedna vrstva
tavg
class       : SpatRaster 
size        : 628, 1068, 1  (nrow, ncol, nlyr)
resolution  : 500, 500  (x, y)
extent      : 263229.1, 797229.1, 5359788, 5673788  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633) 
source(s)   : memory
name        :  sum 
min value   :  1.5 
max value   : 12.4 
# jména vrstev můžeme opravit funkcí terra::names()
names(tavg) <- "tavg"

# podobně existuje funkce terra::time() pro nastavování času u vrstev

# zapišme do GeoTIFF souboru
writeRaster(tavg,
            "outputs/tavg.tif", # máme založenou složku outputs
            overwrite = T) # jinak bychom při opakovaném uložení stejného souboru naráželi na existující soubor

Interakce mezi rastrovými a vektorovými geodaty

Rastry často vystupují v interakcích s vektorovými geodaty, např. se zajímáme jen o nějaké území dané polygonem. Podobně jako u vektorů zde máme funkci crop(). Namísto st_intersection() používáme funkci mask(), nebo rovnou zaměstnáme argument mask ve funkci crop(), pokud je vektorová vrstva, se kterou v obou funkcích pracujeme, totožná.

# řekněme, že se zajímáme jen o budoucí teplotu vzduchu na území Prahy
tavg_praha <- crop(tavg,
                   praha |> 
                     st_transform(crs(tavg)), # musíme si dávat pozor na odlišné crs; současné verze terra by však již měly transformovat vektor na crs rastru automaticky s varováním
                   mask = T)

tm_shape(tavg_praha) + 
  tm_raster(col.scale = tm_scale_continuous(values = "brewer.reds"),
            col.legend = tm_legend(reverse = T))

Velmi častou interakcí mezi rastrem a vektorem je extrahování hodnot buněk rastru podle vektorové vrstvy. U bodů to není nutné, ale u čar a polygonů je vyžadována ještě funkce, kterou chceme hodnoty buněk agregovat. Velmi často je takovou funkcí mean(). Ukažme, jak můžeme extrahovat hodnoty rastru budoucí teploty vzduchu do polygonů jednotlivých krajů Česka.

# vrstvu krajů budeme potřebovat v crs rastru
kraje_t <- kraje |> 
  st_transform(crs(tavg)) |> 
  select(NAZ_CZNUTS3) # ještě raději vybereme podstatný sloupec, aby ostatní nerušily prohlížení výsledných dat

kraje_t <- extract(tavg,
                   kraje_t,
                   fun = mean,
                   bind = T) |> # tento argument zaručí připojení nového sloupce k vektorové vrstvě, což je častý případ, který chceme
  st_as_sf() |> # výsledkem je třída SpatVector, kterou chceme konvertovat na sf
  as_tibble() |> # tady už jen chceme z obyčené data frame dostat tibble
  st_sf() |> # převádíme zpátky na sf
  mutate(tavg = round(tavg, 1)) # a zaokrouhlujeme přidané hodnoty s teplotou

# prohlédneme
kraje_t
Simple feature collection with 14 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 292599.6 ymin: 5377849 xmax: 779122.6 ymax: 5656242
Projected CRS: WGS 84 / UTM zone 33N
# A tibble: 14 × 3
   NAZ_CZNUTS3           tavg                                           geometry
 * <chr>                <dbl>                                      <POLYGON [m]>
 1 Hlavní město Praha    11.3 ((463736.9 5557917, 462832.4 5557917, 462866.1 55…
 2 Středočeský kraj      10.3 ((509972.8 5596562, 510134.3 5596443, 510191.4 55…
 3 Jihočeský kraj         9.2 ((488620.3 5491979, 488509 5492183, 488394.5 5492…
 4 Plzeňský kraj          9.3 ((323887.7 5498025, 323897.2 5497841, 323956.4 54…
 5 Karlovarský kraj       8.5 ((353515.6 5537626, 353506.3 5537641, 353546 5537…
 6 Ústecký kraj           9.8 ((409105.8 5619468, 409098.4 5619421, 409057.5 56…
 7 Liberecký kraj         9   ((500478.1 5639506, 500397.1 5639380, 500211.6 56…
 8 Královéhradecký kraj   9.5 ((570895.7 5612425, 570817.7 5612555, 570696.1 56…
 9 Pardubický kraj        9.7 ((610500.1 5556148, 610477.4 5556166, 610452.5 55…
10 Kraj Vysočina          9.3 ((588390.4 5467659, 588414.7 5467706, 588721.6 54…
11 Jihomoravský kraj     10.9 ((655531.4 5453565, 655707.6 5453753, 655753 5453…
12 Olomoucký kraj         9.7 ((648802.3 5470021, 648771.5 5469937, 648768.3 54…
13 Zlínský kraj          10   ((701371.6 5480437, 701406.5 5480395, 701453.8 54…
14 Moravskoslezský kraj   9.6 ((682666.6 5571834, 682595.4 5571833, 682565.3 55…

Reference

HIJMANS, R. J. (2025): terra: Spatial Data Analysis.
PEBESMA, E. J., BIVAND, R. (2023): Spatial Data Science: With Applications in R. CRC Press, Boca Raton.

Poznámky pod čarou

  1. Občas také proto, že tyto vztahy ignorujeme, uvidíme při aplikaci některých funkcí varování, ale to nám vadit nebude.↩︎

  2. Jde o funkci st_geometry_type(). Dokonce existuje funkce na tázání se na prázdnou geometrii, tj. st_is_empty().↩︎

  3. Aby funkce pracovala správně, potřebujeme mít nainstalovaný balíček ncdf4. Když nebude nainstalovaný, objeví se o tom hláška.↩︎